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Abstract 

The quadratic adaptive integrate-and-fire model Pzhikevich, 2003} 
[izhikevic h, 2007P is recognized as very interesting for its computa- 
tional efficiency and its ability to reproduce many behaviors observed 
in cortical neurons. For this reason it is currently widely used, in par- 
ticular for large scale simulations of neural networks. This model em- 
ulates the dynamics of the membrane potential of a neuron together 
with an adaptation variable. The subthreshold dynamics is governed 
by a two-parameter diff"erential equation, and a spike is emitted when 
the membrane potential variable reaches a given cutoff value. Subse- 
quently the membrane potential is reset, and the adaptation variable 
is added a fixed value called the spike-triggered adaptation parame- 
ter. We show in this note that when the system does not converge to 
an equilibrium point, both variables of the subthreshold dynamical 
system blow up in finite time whatever the parameters of the dynam- 
ics. The cutoff is therefore essential for the model to be well defined 
and simulated. The divergence of the adaptation variable makes the 
system very sensitive to the cutoff: changing this parameter dra- 
matically changes the spike patterns produced. Furthermore from a 
computational viewpoint, the fact that the adaptation variable blows 
up and the very sharp slope it has when the spike is emitted implies 
that the time step of the numerical simulation needs to be very small 
(or adaptive) in order to catch an accurate value of the adaptation 
at the time of the spike. It is not the case for the similar quartic 
( [Touboul, 2008"! ) and exponential ( [Brette and Gerstner, 2005D mod- 
els whose adaptation variable does not blow up in finite time, and 
which are therefore very robust to changes in the cutoff value. 



1 Introduction 

During the past few years, in the neuro-computing community, the problem 
of finding a computationally simple and biologically realistic model of neuron 
has been widely studied, in order to be able to compare experimental record- 
ings with numerical simulations of large-scale brain models. The key problem 
is to find a model of neuron realizing a compromise between its simulation 
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efficiency and its ability to reproduce what is observed at the cell level, of- 



ten considering in- vitro experiments (Koch and Segev, 1998 Izhikevich, 2004 



Rinzel and Ermentrout, 1989). Among the variety of computational neuron 
models, nonlinear spiking models with adaptation have recently been studied 
by several authors ( [Izhikevich, 2004[ [Brette and Gerstner, 20051 |Touboul, 20081 ) 
and seem to stand out. They are relatively simple, i.e. mathematically tractable, 
efficiently implemented, and able to reproduce a large number of electrophysio- 
logical signatures such as bursting or regular spiking. These models satisfy the 
equations: 

ft =Fiv)-w + I 
\^ =a{bv-w) 

where a and h arc non-negative parameters and F{v) is a regular strictly convex 
function satisfying assumption: 

Assumption (Al). There exists e > and a > for which F{v) > av^~^^ 
when V ~> oo (we will say that F grows faster than v^^^ when v ^ (x). 

A spike is emitted at the time t* when the membrane potential v reaches a 
cutoff value 9. At this time, the membrane potential is reset to a constant value 
c and the adaptation variable is updated to w{t*)+d where w{t*) is the value of 
the adaptation variable at the time of the spike and d > is the spike-triggered 
adaptation parameter. 

For these models we prove in section [2] that the membrane potential blows up 
in finite time. Among these models, the quadratic adaptive model Pzhikevich, 2004| 
corresponds to the case where F(v) = , and has been recently used by Eugene 



Izhikevich and coworkers (Izhikevich and Edelman, 2008) in very large scale 



simulations of neural networks. The adaptive exponential model QBrette and Gerstner, 2005| ) 
corresponds to the case where F{v) = e" , has the interest that its parameters can 
be related to electrophysiological quantities, and has been successfully fit to in- 
tracellular recordings of pyramidal cells ( Clopath et al., 20071 jJolivet et al., 2008 ) . 



The quartic model ( jTouboul, 2008P corresponds to the case where F{v) 
v'^ A-2av and has the advantage to of being able to reproduce all the behaviors 
featured by the other two and also self-sustained subthreshold oscillations which 
are of particular interest to model certain nerve cells. 

In these models, the reset mechanism makes critical the value of the adap- 
tation variable at the time of the spike. Indeed, when a spike is emitted at 
time t* , the new initial condition of the system ([T]) is (c, w{t*) + d). Therefore, 
this value governs the subsequent evolution of the membrane potential, and 
hence the spike pattern produced. For instance in ( jTouboul and Brette, 2008al 
jTouboul and Brette, 2008bl ), the authors show that the sequence of reset loca- 
tions after each spike time shapes the spiking signature of the neuron. 

Hence characterizing the reset location of the adaptation variable is essential 
to characterize the spiking properties of these models. To this end, we precisely 
study in this note the orbits of equation ([1]) in the phase plane {v^ w) in order 
to characterize the value of the adaptation variable at the time of the spike. We 
prove in section |2] that the adaptation variable diverges when w — > oo in the 
case of the quadratic model and converges in the cases of the exponential and 
of the quartic model, and study in section [3] the consequences of this fact on the 
spiking signatures and on numerical simulation methods. 
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2 Adaptation variable at the times of the spikes 



As we can see in equation ([1}, the greater the membrane potential the greater 
the derivative of the adaptation variable. When the membrane potential blows 
up, the adaptation variable may either remain bounded or blow up, depending 
on the shape of the divergence of v. When this divergence is not fast enough, 
the adaptation variable simultaneously blows up. 

We prove here that for the models satisfying assumption |(Al)| the membrane 
potential blows up in finite time. We also prove that for quadratic adaptive 
models the adaptation variable blows up at the same time as a logarithmic 
function of v, whereas if there exists e > such that F{v) grows faster than 
■y^+e when w — » cx), then the adaptation variable remains bounded when v —> oo. 

In ( [Touboul, 2008P , we have seen that there exists possibly one stable fixed 
point for system[Tl which corresponds to a resting state. In (ITouboul and Brette, 2008b), 



we prove that all the orbits of the system that do not converge to this stable 
fixed point will be trapped after a finite time in a zone fully included in the half 
space {w < bv} called the spiking zon^. Denote to a time such that the orbit 
is inside the spiking zone. In this zone, we have 

— > F(v) -bv + I 
dt 

It is simple to prove that the solution of the equation 

=F{u) -bu + I 
[u{to) = v{to) 

blows up in finite time under the assumption (Al j^ . Using Gronwall's theorem 
dGronwall, 19191 ) we conclude that v{t) > u{t) and hence v blows up in finite 
time. 

To prove the divergence of the adaptation variable when the membrane 
potential blows up in the case of the quartic model, we study the orbit of 
a solution {v{t)^w{t)) of the differential system ([T]) such that the membrane 
potential blows up at time t* , and characterize the behavior of w{t) in func- 
tion of v{t). In the spiking zone, we have seen that < bv{t) and there- 
fore F{v) — w + I > F{v) ~ bv + I which tends to infinity when v tends 
to infinity. Since v(t) blows up there exists a time ii £ [^0:^*) such that 
we will have F{v{t)) - w{t) + / > fc > for all t G [h.t*). We denote 
(vi := v{ti),'Wi :— w(ti)). After time ti, because of this inequality, the tra- 
jectory in the phase plane can be written as the graph of a function W{v) that 



^We can prove more generally that when F(v)/v'^ tends to a finite constant (possibly 0), 
the adaptation variable will blow up when the membrane potential blows up 

^In the case where the subthreshold system has no fixed point this property can be derived 
from the shape of the vector field in the phase plane, as well as in the case where the initial 
condition (v, w) is such that v is greater than the largest ti- value of the fixed points (the biggest 
solution of F(v) — bv + I = 0) and w < bv: in this case the vector field on the line w = bv 
implies that the trajectory keeps trapped in this zone. In the case where there exist fixed 
points, the proof is slightly more complex and involves the description of the stable manifold 
of the saddle fixed point. 

^For the quadratic model we can get analytic expressions of the solutions involving the 
tangent function, and therefore can derive an upperbound of the explosion time. 
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satisfies the equation: 

■ dW _ g (bv-W) 



{ dW 
dv 



F{v)-W+I 
= Wl 



(i.e. w{t) = W{v{t)) for t G [ti,t*)). Since w{t) is increasing for t G [ti,t*), we 
necessarily have: 

dW ^ ajby-W) 

dv - F{v) -wi+I ^ ' 

Therefore Gronwall's theorem f Gronwall, 1919[ ) ensures us that the solution of 
equation ^ will be lowerbounded for v > vi hy the solution of the linear 
ordinary differential equation: 

( dz^ _ a {b v-z) 
J dv ~ F{v)-wi+I 



Z{vi) = Wl 



that reads: 



\Jv, Fiu) ~wi+I J 

where g{v) — — /J^ F(u)-wi+i - Because of assumption |(A1)[ the integrand 
is integrable, and the function g has a finite limit g{oo) when v — > oo. The 
exponential terms will hence converge when u — > oo. But the integral involved 
in the particular solution diverges in the quadratic case0, since the integrand is 
equivalent when u ^ oo to 

U 

Hence the solution of the linear differential equation @ tends to infinity 
when w — > oo faster than a logarithmic function of v, and so does W{v)^ and 
hence 'w{t) blows up at the time when v{t) blows up. 

Let us now upperbound the adaptation variable on the orbits of the system. 
Using the same notations, since wi < w{t) < bv{t) for t £ [ti,t*), we have: 

dW _^ a{bv~wi) 

dv - F{v)~bv + I ^ ' 

and hence 

X n aibu — Wl) 
W{v) <wi+ -—\ , ' ^ du 



„^ F{u) -bu + I 

In the case where F(u) = v? this integral is bounded by a logarithmic 
function of v and in the case where F{u) grows faster than this integral 

converges when v — > oo. Furthermore, since W is an increasing upperbounded 
function, it converges when w — > oo. 

We therefore conclude that in the case of the quadratic adaptive model, the 
adaptation variable blows up at the explosion time of the membrane poten- 
tial variable v and this divergence is logarithmic in v (see figure [T] (h)), and in 
the case of the quartic and exponential models, the adaptation variable con- 
verges. Note that the adaptation variable diverges whatever the parameters of 



or when F{v) grows slower than v 
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the system provided the membrane potential variable blows up. The value of 
the adaptation variable at the cutoff 9 is simply given by W{6), that depends 
on the parameters of the system and of the initial condition. In the case of the 
quadratic model it is an unbounded increasing function of 9 , and in the quartic 
and exponential models, a converging function of 9. 



3 Consequences 

The divergence of the adaptation variable at the times of the spikes significantly 
impacts the theoretical, qualitative and computational analysis of the model. It 
appears to be a critical parameter of the quadratic model. 

We have seen that changing the cutoff value resulted in changing the value 
of the adaptation variable at the times of the spikes. Let (wo,wo) be an initial 
condition for the system ([T]). If the neuron fires, its membrane potential will 
reach the cutoff value at a given time. Since the membrane potential blows 
up in finite time, the time of the first spike emitted will not be very sensitive 
to changes in the cutoff value provided it is high enough. But the after-spike 
reset location (c, W{9) +d) will significantly change when varying 9. The whole 
subsequent evolution of the system is therefore affected, as soon as the second 
spike is emitted. Thus the spike pattern produced depends on the cutoff value. 

In the case of the quartic and exponential models, the adaptation variable 
converges when the cutoff tends to infinity. Therefore, the model defined by ([T]) 
with an infinite cutoff value is mathematically well defined. In that spike 
is emitted when the membrane potential blows up and subsequently we reset the 
membrane potential to a fixed value c and add to the value of adaptation vari- 
able at the explosion time the spike-triggered adaptation parameter. We call this 
system the intrinsic system. The behavior of the system and the spike patterns 



it produces can be mathematically studied (see (Touboul and Brette, 2008a 



Touboul and Brette, 2008bl )). Interestingly, these intrinsic spike patterns un- 



dergo bifurcations with respect to the parameters of the model. When consid- 
ering a finite cutoff, the model (or the numerical simulation) will approximate 
these intrinsic behaviors provided that the cutoff threshold is high enough. The 
sensitivity to the cutoff in these cases will hence be very limited except in very 
small regions of the parameter space around the bifurcations of the intrinsic 
system. Unfortunately, for the quadratic model, no intrinsic behavior can be 
defined because of the divergence of the adaptation variable: the behaviors it 
produces will depend on the choice of the threshold. 

First of all, we have seen that the dependency of this reset location in the 
quadratic model is a logarithmic function of 0, which makes the variations of the 
reset value in function of the cutoff unbounded but quite slow. Small changes 
in the cutoff slightly impact the value of the reset adaptation variable. For 
instance if we consider the firing rate of a neuron in the case where the system 
has no fixed point, increasing the cutoff value results in the case of the quadratic 
model in a slow continuous decrease of the firing rate of the neuron that tends 
to zero as the cutoff increases, whereas the firing rate converges for the quartic 
model to the related intrinsic firing rate (see figure [TJ(g)). 

When considering the spike patterns produced, the effects of changes in 
the cutoff value for the quadratic model are much more dramatic. Indeed, 
the sequence of adaptation values at the times of spikes shapes the spike pat- 
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tern produced: for instance, regular spiking corresponds to the convergence 
of this sequence, and bursting to cycles in this sequence. These properties 
are very sensitive to changes in the parameters of the model: bifurcations 
between different spike patterns, and even chaos appear when the model's 
parameters vary (see ( |Touboul and Brette, 2008a| [Touboul and Brette, 2008bl 
|Naud et al., 2008| ). In the case of the quadratic model, we have seen that these 
adaptation values strongly depend on the cutoff. Therefore, since the depen- 
dency on the cutoff is unbounded, from a given initial condition and for fixed 
values of the parameters, increasing the cutoff may result in crossing many bi- 
furcation lines, and hence in producing many different behaviors. We present 
in figure [1] a graph showing that bifurcations and chaos occur with respect to 
the cutoff value, in the usual range of simulation parameters. For instance, a 
period doubling bifurcation appears when varying the cutoff value (in figure[Ue) 
we give a graph of the stationary reset values in function of the threshold 9), 
that results in abruptly switching from a regular spiking behavior to a bursting 
behavior (figures [TJa) and[ljb)). More complex bifurcation structures involving 
chaotic patterns also appear, and in this case, infinitesimal changes in the cutoff 
value result in dramatic changes in the behavior. This raises the question of 
the meaning of the cutoff value in these ranges of parameters (see figure [TJf)). 
Changing the cutoff in that case makes the system switch between chaotic spik- 
ing, bursts with 8, 4 and eventually 2 spikes, for the cutoff values considered. 
And this behavior will not be observed only for very particular values of the 
parameters of the system. Depending on the extension of the interval where the 
cutoff value varies, quite a large set of parameters will present bifurcations in 
the nature of the emitted spike train. 

Because of this sensitivity, the cutoff value and the different parameters of the 
model have to be very carefully evaluated in order to quantitatively fit datasets. 
In this context the meaning of the threshold and therefore the problem of its 
accurate evaluation has to be specifically addressed in the case of the quadratic 
model, since it has no clear biophysical interpretation. 

Eventually, from the numerical viewpoint, the unboundedness of the adap- 
tation variable and of its time derivative at the explosion times of the membrane 
potential makes the accurate computation of this value very difficult. In par- 
ticular, the time step necessary to accurately estimate this value has to be very 
small (or to be adaptive as a function of the value of the membrane potential 
variable) in order to obtain the right spike pattern. These remarks relativize 
the statement that this model can be efficiently simulated since very accurate 
methods have to be implemented in order to correctly evaluate the adaptation 
variable at the time of the spike. 

These remarks do not apply for the models where the adaptation variable 
converges at the times of the spikes. In these cases, the system has intrinsic 
properties that make the times of the spike and the adaptation variable at these 
times robust to the choice of the cutoff value provided it is big enough and the 
numerical simulations less sensitive to the choice of the time step. However, 
the exponential integrate-and-fire model involves an exponential function which 
increases very fast as a function of the membrane potential. This fact implies 
that the time step chosen for the simulation has to be very small and makes 
this system less suitable for large scale simulations. In particular, the numerical 
scheme must be designed carefully in order to capture the value of the adaptation 
variable at the times of the spike. Simulating the orbit in the phase plane for 
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Figure 1: Sensitivity of the spike patterns with respect to the cutoff value for 
the quadratic model, for different set of parameters. Parameters used: (A) = 
{a = 0.02; & = 0.19; c = -60; d = 1.419,7 10.25}; (B) = {a = 0.1; 6 = 
0.26; c = -60; d = 0;}, (C) = {a = 0.02,6 = 0.19, c = -57.7, d = 1.15,/ = 
10.377}. Simulations done with an Euler scheme with time steps ranging from 
10"'* to 10~^. For figure (a) and (b) the parameters used arc (A) with cutoff 
of 30 and 45 respectively: a small increase of the cutoff results in a sharp 
transition from spiking to bursting linked with a period doubling bifurcation 
for the adaptation value at the reset represented in figure (e). Since there is a 
bifurcation, the qualitative change is very sharp. Figures (c) and (d) corresponds 
to the parameters (B) with cutoffs value 32.9 and 33 respectively. Changing the 
cutoff results in two very different global behaviors. Fig. (c) and (f) represent 
the stationary sequence of reset values as functions of the threshold 6. Figure (f) 
corresponds to the set of parameters (C) for cutoff values ranging from 20 to 100: 
an intricate bifurcation structure appears. Figure (g) shows the convergence of 
the firing rate to the intrinsic firing rate in the case of the quartic model, while 
the firing rate of the quadratic model regularly decreases to 0. (h): Divergence 
of the variable w in function of v in semilogarithm axis. 
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instance would result in computing very accurately this essential value. 

Conclusion 

In this note we proved that the adaptation variable of the adaptive quadratic 

model blew up at the times of the spikes whereas it converged for the quartic 
and the adaptive exponential models. This property has some important im- 
plications that are discussed in the note. Prom a theoretical point of view, we 
showed that the nature of the spike patterns produced undergoes bifurcations 
with respect to the cutoff value, and this made the system very sensitive to 
this parameter: small changes in the value of this parameter can deeply affect 
the nature of the spiking pattern. From a quantitative viewpoint, it raises the 
question of how to evaluate this threshold in order to fit datasets, and from a 
numerical viewpoint, it has implications on the efficiency of the simulation al- 
gorithms to use. The convergence of this value for models having a faster blow 
up at the times of the spike, such as the quartic or the exponential adaptive 
models, implies that the system presents intrinsic spiking properties which can 
be mathematically studied, efficiently simulated and robust to changes in the 
cutoff value. 

Acknowledgments 

The author warmly acknowledges Bruno Cessac, Olivier Faugeras and Frangois 
Grimbert for interesting discussions and suggestions. 

References 

Brette, R. and Gerstner, W. (2005). Adaptive exponential integrate-and-fire 
model as an effective description of neuronal activity. Journal of Neuro- 
physiology, 94:3637 3642. 

Clopath, C, Jolivet, R., Ranch, A., Liischer, H., and Gerstner, W. (2007). 
Predicting neuronal activity with simple models of the threshold type: 
Adaptive Exponential Integrate-and-Fire model with two compartments. 
Neurocomputing, 70(10-12):1668-1673. 

Gronwall, T. H. (1919). Note on the derivative with respect to a parameter of 
the solutions of a system of differential equations. Annals of Mathematics^ 
20:292-296. 

Izhikevich, E. (2003). Simple model of spiking neurons. IEEE Transactions on 
Neural Networks, 14(6): 1569-1572. 

Izhikevich, E. (2004) . Which model to use for cortical spiking neurons? IEEE 
Trans Neural Netw, 15(5):1063-1070. 

Izhikevich, E. (2007). Dynamical Systems in Neuroscience: The Geometry of 
Excitability And Bursting. MIT Press. 



8 



Izhikevich, E. M. and Edelman, G. M. (2008). Large-scale model of mammalian 
thalamocortical systems. Proc Natl Acad Sci USA, 105(9) :3593-3598. 



Jolivet, R., Kobayashi, R., Rauch. A., Naud, R., Shinomoto, S., and Gcrstner, 
W. (2008). A benchmark test for a quantitative assessment of simple neuron 
models. Journal of Neuroscience Methods, 169(2) :417-424. 

Koch, C. and Segev, I., editors (1998). Methods in Neuronal Modeling: From 
Ions to Networks. The MIT Press. 

Naud, R., Macille, N., Clopath, C, and Gerstner, W. (2008). Firing patterns in 
the adaptive exponential integrate-and-fire model. Biological Cybernetics 

(submitted). 

Rinzel, J. and Ermentrout, B. (1989). Analysis of neural excitability and oscil- 
lations. MIT Press. 

Touboul, J. (2008). Bifurcation analysis of a general class of nonlinear integrate- 
and-fire neurons. SIAM Journal on Applied Mathematics, 68(4):1045-1079. 

Touboul, J. and Brette, R. (2008a). Dynamics and bifurcations of the adaptive 
exponential integrate-and-fire model. Biological Cybernetics (submitted). 

Touboul, J. and Brcttc, R. (2008b). Spiking dynamics of bidimensional 
integrate-and-fire neurons. Research Report 6531, INRIA. 



9 



